Vacancy induced magnetism in graphene and graphene ribbons 
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We address the electronic structure and magnetic properties of vacancies and voids both in 
graphene and graphene ribbons. Using a mean field Hubbard model, we study the appearance 
of magnetic textures associated to removing a single atom (vacancy) and multiple adjacent atoms 
(voids) as well as the magnetic interactions between them. A simple set of rules, based upon Lieb 
theorem, link the atomic structure and the spatial arrangement of the defects to the emerging 
magnetic order. The total spin S of a given defect depends on its sublattice imbalance, but some 
defects with S — can still have local magnetic moments. The sublattice imbalance also determines 
whether the defects interact ferromagnetically or antiferromagnetically with one another and the 
range of these magnetic interactions is studied in some simple cases. We find that in semiconduct- 
ing armchair ribbons and two-dimensional graphene without global sublattice imbalance there is 
maximum defect density above which local magnetization disappears. Interestingly, the electronic 
properties of semiconducting graphene ribbons with uncoupled local moments are very similar to 
those of diluted magnetic semiconductors, presenting giant Zeeman splitting. 

PACS numbers: 



I. INTRODUCTION 

Magnetic order occurs, in most instances, in materials 
with partially filled d or / shells. There is, however, a re- 
cent awareness that the possibility of magnetic order can 
also occur in materials without open d or / shells 1 -^! 4 -! 5 -. 
Experimental evidence of this new type of magnetism has 
been found in thin films of certain oxides (Hf02, ZnO, 
TiC^) 1 , as well as irradiated graphite 2 , and thiol-capped 
gold nanoparticles^ 7 -. 

Although more experimental work is probably neces- 
sary to confirm and understand magnetism in these sys- 
tems, there are at least two scenarios for which theory 
provides a mechanism for the appearence of magnetism 
without d or / open shells. On one side, in some lat- 
tices intrinsic lattice defects like vacancies lead to the 
formation local magnetic moments, a preliminary condi- 
tion for the existence of magnetic order. This is the case 
in graphite 8 -, graphene 9 - ' 10 ' 11 and II- VI semiconductors 1 ^. 
On the other side, it has been recently found that clusters 
with specific shapes, like triangular graphene islands 3 - or 
icosahedral 4 ' 5 gold clusters, have large degeneracies at 
the Fermi energy in their single particle spectra. These 
degeneracies are related to the symmetry of the nanos- 
tructure and, in words of Luo et a/£, they behave like 
'superatoms', with magnetic ground states that comply 
with atomic-like Hund's rules. 

Importantly, both vacancy induce d 10 ! 11 and 'super- 
atomic' magnetism 5 occur in graphene structures and, 
as we show here, have the same origin. In this work 
we present extensive numerical work to understand va- 
cancy induced magnetism in graphene and graphene rib- 
bons and we analyze our results in the context of a 



broader theoretical framework that unifies superatomic 3 
and vacancy induced magnetis m 10 ! 11 in graphene. In part 
our motivation stems from the recently shown possibility 
of fabricating high mobility graphene based field effect 
transistor s 13 ' 14 ! 15 ' 16 ! 17 ' 18 which has created enormous in- 
terest in graphene based electronics^ 9 .. Additional pos- 
sibilities arise from the fabrication of semiconducting 
graphene ribbon s 20 ' 21 ' 22 and graphene nanoisland o 19 ' 23 
with top-down techniques as well as the growth of 
graphene islands with bottom-up technique s 24 ' 25 . Mag- 
netic order in patterned or nanostructured graphene 
would bring up new opportunities of research in spin- 
tronics. 

Graphene honeycomb structure is a bipartite lattice, 
formed by two interpenetrating triangular sublattices, A 
and B, such that the first neighbors of an atom A belong 
to the sublattice B and viceversa^ 6 -. The low-energy elec- 
tronic structure of graphene can be described by a single- 
orbital (p z ) nearest-neighbor hopping Hamiltonia n 26 ' 27 . 
This model correctly describes two dimensional graphene 
as a zero gap semiconductor with linear bands around the 
Fermi energy. The single particle spectrum of a nearest- 
neighbor tight-binding model in a bipartite lattice has 
particle-hole symmetr y 28 ' 29 . 

The magnetic properties of both graphene-based 
nanostructures and defective graphene are intimately re- 
lated to the appearance of midgap states and how they 
are affected by electron-electron interactions. The exis- 
tence of zero-energy states in disordered bipartite lattices 
was proved by Inui et al. 28 . Within the first-neighbor 
tight-binding model, a sufficient condition 3 -^ for the ex- 
istence of midgap states is the existence of a finite sub- 
lattice imbalance, Nj = Na — Nb, where Na and Nb are 
the number of atoms belonging to each sublattice or miss- 
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ing from each sublattice in an otherwise perfect system. 
Thus, whereas ideal graphene has Nj — and no midgap 
states, both defective graphene and some graphene is- 
lands, such as triangles, can present finite sublattice im- 
balance and |iV/| midgap states. The result of Inui et al3& 
has been used in a recent work on disorder in graphene by 
Pereira et al£&. Incidentally, the existence of zero-energy 
or midgap states in uncompensated graphene structures 
was known long ago in the context of chemical studies of 
hydrocarbons as the Longuet-Higgis conjecture 3 -^. 

Because of the particle-hole symmetry, midgap states 
are half filled for neutral graphene and the appearence of 
magnetic moments is expected in analogy with Hund's 
rule in atomic magnetism. The Hubbard model extends 
the single-particle tight-binding model including the ef- 
fect of Coulomb repulsion between two electrons in the 
same atomic site. Importantly, a theorem by Lieb, valid 
for the exact ground state of the Hubbard model and 
neutral bipartite lattices 3 -^ states that the total spin S of 
the ground state is given by 2S — \Na — Nb\ = \Ni\. 
Lieb's theorem provides a rigorous connection between 
vacancies in the graphene lattice and the emergence of 
magnetism. As as result, sublattice unbalanced neutral 
graphene will always present a finite total magnetic mo- 
ment. 

Although Lieb's theorem provides the total spin of the 
ground state, it does not say much about the actual lo- 
cal magnetic order or spin texture. For instance, S = 
does not preclude the existence of local magnetic mo- 
ments coupled antiferromagnetically or presenting com- 
pensated ferrimagnetic order. The most notorious ex- 
ample of compensated ferrimagnetic order can be found 
in zigzag ribbon s 33 ' 34 ! 35 ' 36 ' 37 , where each edge presents 
ferromagnetic order antiparallel to each other for a to- 
tal vanishing magnetic moment. Other examples can be 
found in hexagonal graphene islands, where, beyond a 
critical size, contiguous sides alternate the direction of 
the ferromagnetically ordered magnetic moments 3 . 

The rest of this paper is organized as follows. In Sec. 
im we review the single orbital Hubbard model and the 
different methodologies used to describe the electronic 
structure of defective graphene and graphene ribbons. 
The underlying non-interacting spectrum and associated 
magnetic textures can be anticipated following some ba- 
sic rules which are presented in Sec. IIIIl We illustrate 
the validity of the rules by numerical calculations in the 
case of semiconducting armchair ribbons (Sec. IIV|) with 
vacancies, voids, or notches, both in the non-interacting 
(Sec. IIV[) and interacting (JVj cases. The results for bulk 
graphene are discussed in Sec. (|VTJ). Summary and con- 
clusions are presented in Sec. I VIII 



II. METHODOLOGY 

We consider the low-energy physics that takes place 
in the subspace expanded only by the single p z orbital 
(the one perpendicular to the graphene plane). Next- 



to-near neighbor hopping is neglected and the electron- 
electron interactions are included locally in the form of 
an on-site repulsion or Hubbard model. When the inter- 
actions are turned off this reduces to the widespread one 
orbital tight-binding modelSli 3 ^^. The Hubbard term 
is treated in a mean field approximation^ 3 -. Comparison 
between the results so obtained and density functional 
theory (DFT) calculations yield very good agreement 
for two-dimensional graphene 4 ^, carbon nanotube o 41 ' 42 , 
zigza g 35 ! 37 and armchair graphene ribbons^ 4 -, as well as 
graphene islands 3 -. 

We model vacancies and voids in perfect graphene or 
graphene ribbons by removing atoms, actually, by re- 
moving the representing p z orbitals in the tight-binding 
model. This results in a reduction of the coordination of 
the atoms adjacent to the missing atoms. We ignore the 
lattice distorsion and we assume that the on-site energy 
is the same for edge and bulk atoms. The single-orbital 
hamiltonian implicitly assumes full hydrogen passivation 
of the sp2 dangling bonds of the atoms without full co- 
ordination. This assumption, which might not be com- 
pletely realistic in the case of actual vacancies^, does not 
invalidate our model; for we can consider an alternative 
physical realization: The chemisorption of a hydrogen 
atom on top of a bulk graphene atomii effectively re- 
moves a p z orbital from the low-energy hamiltonian. In 
our one-orbital model there is no difference between these 
two scenarios. 

DFT calculations on graphene ribbon s 34 ! 36 , graphene 
island s 3 ' 46 ! 50 , and bulk graphene with vacancies^ 1 - have 
shown that the results follow the predictions of Lieb's 
theorem, even though DFT calculations go beyond the 
first-neighbor hopping, short-range interaction Hubbard 
model on which the theorem is based. In other words, 
second neighbour hopping and inter-site Coulomb repul- 
sion, present in the DFT calculations, do not modify the 
relation between lattice imbalance and total spin of the 
ground state warranted for the Hubbard model for which 
these couplings are absent. From this point of view these 
couplings are irrelevant. It is thus justified to consider 
the following mean-field hamiltonian: 

H = H Q + U ^2{nii(nu) + ni;(n iT )) - U ^(?i u )(n lT ), 

i i 

(i) 

where i runs over all lattice sites and the non-interacting 
hamiltonian reads 

ffo = £>(ctc 7 -+cjc i ), (2) 

where the sum runs over nearest- neiborgh lattice sites i, j 
and t = 2.5 eV. Without loss of generality, we have set the 
diagonal terms of the Hamiltoninan to zero. For neutral 
graphene we can rewrite the mean-field hamiltonian (up 
to a constant) as the sum of two terms: 

H = H a + ^J2 n *( n <) -U^22m l {m l ) (3) 

i i 
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where mt = \ (n,-f — n,jj and m — + mi- The sec- 
ond term in H represents the non-trivial contribution of 
interactions. 

The calculations have been performed considering 
three different types of boundary conditions. For the 
evaluation of the non-interacting density of states (DOS) 
in ribbons we compute the Green's function projected on 
the region where the defects are located. The perfect re- 
gions of the ribbon away from the defects are included in 
the Green's function by means of a self-energy. We refer 
the reader to Ref. 40; for more details on this method- 
ology. When the interactions are turned on (U ^ 0) we 
consider ribbons with periodic boundary conditions in 
one direction or, for bulk graphene, periodic boundary 
conditions in both directions. More details will be given 
in the respective sections. 

III. BASIC RULES 

In this section we provide some general rules to un- 
derstand the appearance of midgap states and magnetic 
textures due to single-atom vancacies or, more gener- 
ally, voids in an otherwise sublattice balanced graphene 
structure, for instance, an infinite defect-free semicon- 
ducting armchair ribbon. A void in graphene can be 
partially characterized by the number and type, A or B, 
of atoms removed from the otherwise perfect structure. 
We will label voids as one would do for chemical com- 
pounds, Ajv 4 Bjv B . Voids will be unbalanced when the 
are created by removing Na and Nb atoms such that 
Nj = Na — Nb ^ The sublattice imbalance Nj, which 
can be either positive or negative, can be interpreted as 
an imbalance "charge". This quantity is central to our 
discussion, although the exact formula of the void is also 
important since it gives an idea of the size and shape of 
the void. In the case of ribbons the voids can be close 
to the edges, becoming thus notches. (For the most part 
we will refer to voids in the bulk of the ribbon, but the 
conclusions apply equally to the case of notches on the 
edges.) For a single unbalanced void characterized by 
Ni, \Ni\ zero-energy states appear in the non- interacting 
spectrum with weight on only one sublattice^. When the 
graphene structure presents a gap E g , as in armchair rib- 
bons, these states are normalizable and localized around 
the void (in contrast to two-dimensional graphene where 
the zero-energy states are not normalizable 51 ). Figure [T] 
shows various examples of voids in a ribbon with different 
sublattice imbalance Nj. 

Let us consider now two voids, characterized locally by 
Ni(l) and Ni(2), sufficiently separated so that they do 
not affect each other. According to the result of Inui et 
ali2&, the single particle spectrum has, at least, 

Nf a = \Nj(l) + N X (2)\ (4) 

midgap states. The important question is what happens 
when the distance between them is such that they do 
affect each other. If Nj(l) and Nj(2) have the same 




FIG. 1: Examples of voids with different sublattice imbalance 
in the middle of a graphene ribbon. From left to right the 
associated imbalance charge is Ni =0,-1 (vacancy), and 2. 

sign, there are |iVr(l)| + 1^/(2)1 midgap states, regard- 
less of the distance. If they have different signs, e.g., 
Nj(l) + Nj(2) = 0, Eq. 0] apparently warrants the an- 
nihilation of midgap states. Within the non-interacting 
model midgap states are 100 percent sublattice polar- 
ized. The non-interacting Hamiltonian has finite matrix 
elements between states that have weight on different 
sublattices. Hence, the mechanism for midgap state anni- 
hilation is hybridization of midgap states localized in dif- 
ferent sublattices. This annihilation occurs as bonding- 
antibonding pairs of midgap states form, resulting in a 
shift in their energy and in a loss of the sublattice polar- 
ization. For large distances, however, this annihilation 
does not occur. 

A well understood related example occurs in zigzag 
ribbon a 38 ' 39 . The edge of a zigzag ribbon has a local 
sublattice imbalance. If the top edge belongs to the A 
sublattice, the bottom edge belongs to the B sublattice. 
States fully localized in the edge have zero energy and 
are localized in a single sublattice. To the extent that 
states mostly localized on the top edge penetrate into 
the ribbon, they hybridize with states mostly localized 
on the bottom edge. This mixing results in a bonding- 
antibonding splitting that takes these states away from 
the Dirac energy. In the case of zigzag ribbons, the degree 
of localization in the unit cell depends on the wavevector. 
States close to the zone boundary are very localized in 
the edges and have energy very close to zero 38 ! 39 . The 
localization decreases as the wavevector departs from the 
boundary, resulting in the hybridization and the depar- 
ture from zero energy 3 - 8 -. 

In the general case one can conclude that the minimum 
number of zero-energy states will be given by 

Nr n = Y / \ Nl (a)+N I (P)\, (5) 

a,/3 

where the integer indeces a and f3 run over voids with the 
same imbalance sign, respectively. In practice, within an 
arbitrarily small energy interval \E\ — > 0, the number of 
zero-energy states can be as large as 

a (3 

In the general case, Nz will be a number between TV™" 
and Nf ax . 
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When electron-electron interactions are turned on, at 
least locally in the form of a Hubbard-type interaction, 
Lieb's theorem guarantees that |JV/| = 25 for neutral 
graphene. The theorem, however, does not exclude the 
possibility of spin-symmetry broken local magnetic order 
when 5 = or small. For instance, two or more voids 
with local sublattice imbalances that cancel out the total 
imbalance can still retain their local magnetic order when 
they are not in proximity. When the imbalance of the 
void is zero but the size is large an internal ferrimagnetic 
order cannot be discarded either. In general, calculations 
will be necessary to ascertain the spin texture in these 
situations. A few conclusions, however, can be reached 
without actually performing any calculations. One can 
distinguish four cases: 

1. Ng in = N z = N z lax : In this case all the voids 
are of the same sign. The coupling between them 
is always ferromagnetic and the spin of the ground 
state is 25* = Nz- The splitting with smaller spin 
states will depend on the inter-void coupling. 

2. N™ in = N z < N™ ax : In this case all the voids of 
different sign are in proximity and interact, yielding 
a 25 — N™ n state. Calculations will be necessary 
to ascertain the spin texture in these situations. 

3. Nf in < N z = N% ax : In this case all the voids of 
different type are separated and uncoupled. The 
ground state has 25 = A™", but the spin- flip gap 
is negligible since there are uncoupled magnetic mo- 
ments. 

4. N^ lin < N z < N™ ax : In this case there are voids 
of different sign, but some of them are uncoupled 
and some not. This is the most general case. The 
ground state has 25 = A™", but, as in the previ- 
ous case, the spin-flip gap is negligible since there 
are uncoupled magnetic moments. Calculations 
will be necessary to ascertain the spin texture in 
these situations as well. 



IV. DEFECTS IN SEMICONDUCTING 
GRAPHENE RIBBONS: NON-INTERACTING 
THEORY 

In this section we study the electronic structure of de- 
fective graphene ribbons within the non-interacting tight- 
binding model. The results are obtained using the cluster 
embedded method described in Ref. In this method- 
ology a finite portion of the ribbon containing the de- 
fects is attached to two semi-infinite perfect ribbons of 
the same width and compute the DOS of the defective 
region. In the defect free case we obtain a gap in the 
DOS. We consider armchair graphene ribbons of width 
W — N y a, where N y is an integer number and a — 2.42 
Ais the graphene lattice parameter. We only consider 
values of W such that, within the first-neighbor tight- 
binding model, the ribbon is semiconducting 38,39 . This 



happens if N y + 1 is not a multiple of three. More realis- 
tic calculations— predict that, because of lattice distor- 
sion of the edge atoms, even ribbons with N y + 1 = 3m, 
where m is an integer, are semiconductors. Semiconduct- 
ing graphene ribbons attract interest due to possible ap- 
plications in nanoelectronic a 20 i 21 ' 43 i 44 . As in the case of 
Si based semiconductors, their electronic structure might 
be strongly influenced by impurities. Here we study the 
effect of vacancies and voids, which are expected to act 
as neutral impurities. 



A. Single void 

The simplest defective structure is a perfect semicon- 
ducting graphene ribbon from which a single atom, A or 
B, is removed. In agreement to Eq. |4l a zero energy state 
appears in the DOS. For neutral graphene, this state is 
half filled. In other words, an spin unpaired electron oc- 
cupies the midgap state. The local density of states at 
zero energy, which is nothing but the modulus square of 
the wave function associated with the zero energy state, 
is shown in the inset of Fig. [Ha) . The state is localized 
in the neighborhood of the vacancy. The shape of the 
midgap state is also peculiar: It has a clear directional- 
ity. Monoatomic vacancies have the shape of a triangle. 
Two vertices of the triangle point towards the edges of 
the ribbon, whereas the lateral vertex can point down- 
stream or upstream along the ribbon. The midgap state 
is peaked around the lateral vertex. Hence, midgap states 
have a strong directional character. 

Importantly, the integrated charge, including both mid- 
gap and band states below the Fermi energy, yields a ho- 
mogeneous charge distribution: there is one electron per 
atom in every atom, even in the presence of the vacancy. 
Hence, the localized midgap state does not imply charge 
localization, yet there is a finite spin density. The total 
spin 5 of the neutral graphene with one vacancy is 1/2 
and the spin density does show a non-homogeneous tex- 
ture, as shown in Fig. HJa). Hence, the region of the 
material around the defect has no local charge but has 
local spin. 

Our next step is to consider larger voids. In Fig. (2^b) 
we show the results for a void with N] — 2, obtained 
by removing four atoms (A3B1). In agreement to Eq. 
HI there are two midgap states (per spin). Their local 
density of states is shown in the inset of Fig. [2jb) . As 
in the case of monoatomic vacancy, two electrons occupy 
the midgap states. The integrated local charge is also 
homogeneous: One electron per site. Within the frame- 
work of the non-interacting model we cannot discrimi- 
nate between the 5 = or the 5=1. As we discuss 
below, when Hubbard repulsion is turned on Lieb's theo- 
rem warrants that the spin of this structure is 5 = 1. In 
Fig. [UJb) we show the magnetization density, calculated 
including the interactions, as discussed below. As in the 
case of a sigle missing atom, there is a magnetic texture 
with 5=1 which is localized in a region without local- 
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FIG. 2: (Color online) (a) Magnetic moments on lattice sites 
around a single vacancy. Inset: Probability density of the 
zero-enery state built with gaussian functions on lattice sites, 
(b) Same as in (a) but for a triangular void with Ni = 2. 
Inset: Same as in (b) but summing over the two zero-energy 
states. 

ization of extra charge. Triangular voids maximize \Nj\ 
while removing the minimum number of atoms. Larger 
ones (|-/Vr| > 2) exhibit similar features to the ones al- 
ready discussed. More complicated voids with zig-zag 
edges such as hexagons or rombhi, which have Ni = 0, 
can still present quasi zero-energy states if they are suf- 
ficiently big and therefore might exhibit spin textures as 
discussed below. 



B. Two voids 

As a step towards understanding the electronic struc- 
ture of graphene with a finite density of defects, we first 
consider the electronic structure of two voids with the 
same absolute value of A/. Each void has a well defined 
sublattice imbalance number Ni when apart, which can 
be positive or negative. If the sublattice imbalance of 
the two voids has the same sign, the global structure has 
twice as many zero energy states as the separated defects. 
The non-interacting hamiltonian does not couple sites on 
the same sub-lattices so that the zero-energy states as- 
sociated with two vacancies with sublattice imbalance of 
the same sign cannot interact, regardless of the distance 
separating them, i.e., the non-interacting DOS will al- 
ways present a delta function at zero energy that can 
accomodate 2 x \Ni\ electrons per spin channel. 

When the imbalance numbers are of different sign they 
cancel out each other. When the defects are far away 
from each other their local electronic structure is ex- 
pected to be the same as that of a single defect: Midgap 
states localized in a single sublattice around the missing 
atoms. As the defects become closer, the single parti- 



cle Hamiltonian, which couples atoms of different sub- 
lattices, will hybridize the otherwise sublattice polarized 
midgap states, which will result in bonding and antibond- 
ing pairs away from zero energy. The localization length 
of the single defect states sets the length scale at which 
this hybridization occurs. 
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FIG. 3: (Color online) Density of states near the Dirac point 
for an armchair ribbon of W = 7a with two vacancies pre- 
senting charge of different sign and same modulus (Ni = ±1). 
Solid lines correspond to the B+A case (left lower inset) and 
dashed lines correspond to the A+B case (left upper inset). 
A finite broadening has been added for visibility's sake of 
the delta functions. The finite, but small, energy splitting 
in the former case is not visible for this broadening. Right 
inset: Bonding-antibonding energy splitting as a function of 
the distance between vacancies for the two different spatial 
orderings. 

Our numerical calculations confirm this scenario. For 
a given width, the hybridization depends on the distance 
and, given the directional character of the midgap states 
in ribbons, on the relative orientation. In Fig. [3]we show 
the DOS for a system with two monoatomic vacancies A 
and B (Nj = ±1, respectively). They are aligned along 
the ribbon axis and placed at a distance of 6.35a away 
from each other for the two possible spatial orderings, 
A+B (head to head) and B+A (tail to tail) as shown 
in the left insets. Due to the high directional charac- 
ter of the associated zero-energy states, the coupling is 
not invariant against the interchange of positions and the 
zero-energy states hybridize differently, depending on the 
spatial ordering. In one case the two-fold zero-energy 
peak clearly splits into two above and below the Fermi 
energy. In the other the splitting is much smaller (not 
visible in this scale). For one relative orientation the 
wave functions overlap and the degeneracy is strongly 
removed. For the other the wave functions do not couple 
at this distance and the degeneracy is practically unaf- 
fected. In the right inset of Fig. [3]we show a logarithmic 
plot of the energy splitting as a function of the distance 
for the two cases. The splitting decays exponentially in 
both, reflecting the localized character of the zero-energy 
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states. 

We now consider the case of pairs of defects with larger 
sublattice imbalance. Figure 2] shows the DOS for two 
triangular voids characterized by Ni = +2 and Nj = — 2 
(A3B1 and A1B3, respectively) at different distances. We 
have selected only one possible ordering in this case (tail 
to tail). According to these sublattice imbalances each 
void has associated two localized states. These two states 
also present a strong directional character, but is differ- 
ent for the two. This can be inferred from the two dif- 
ferent bonding-antiboding splitting energies for a given 
distance seen in Fig. H We note that, even when the 
voids approach each other, the splitting associated with 
one of the localized states remains small, still being prac- 
tically zero for small distances Only in the extreme limit 
of zero distance when the two voids merge into a single 
one with sublattice imbalance Nj — (A4B4) there are 
no zero-energy states. 
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the ribbon, we compute the non-interacting DOS for an 
A6B4 void plus an AB2 void placed on the edges (i.e., 
notches) with Nj = 2 and Nj — — f, respectively. Re- 
moving just one atom to create a notch with Ni = — 1 
would have given the same charge as the AB2 defect, but 
it would be chemically very unstable and we ignore that 
possibility. The notches are located on opposite edges, al- 
though the results apply the same for notches on the same 
edge. A single doubly-degenerate state appears at zero 
energy for the Ni = 2 notch. When the second notch is 
added in close proximity only a single zero-energy state 
remains, according to the total sublattice imbalance of 
the system Nj = 2 — f = f . 
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FIG. 5: (Color online) DOS for a single notch with imbalance 
charge Ni = 2 (blue dashed line, right inset). The same notch 
with an additional notch nearby of charge Ni = — 1 (black 
solid line, left inset). 



FIG. 4: (Color online) DOS projected on the vicinity of two 
triangular voids placed along the axis of a semiconducting 
ribbon for different relative distances. Fhe imbalance charges 
are the same, but differ in sign (Ni = ±2). 

As the sublattice imbalance of the merging voids be- 
comes bigger and these condense into even bigger Nj = 
voids, the number of states that appear in a vicinity of 
zero \E\ — > increases with the charge of these. Since 
the appearance of magnetic order relies on the existence 
of zero-energy states, large voids with Ni = can still 
present ferrimagnetic order, the only condition being that 
they are formed out of voids with large sublattice imbal- 
ance. In other words, their contours must present suf- 
ficiently long zig-zag sections. This limits the possible 
shapes of these voids to, e.g., rombohedral (see Fig. [5]) 
or hexagonal forms. This conclusion is no different from 
that reached on graphene hexagonal islands 3 or finite 
length ribbons^, where calculations have revealed com- 
pensated ferrimagnetic order developing along the edge 
beyond a certain critical size. 

Finally, in order to stress the fact that there is nothing 
in the previous discussion specific to voids in the bulk of 



In summary, defective structures with sublattice imbal- 
ance result in half-filled midgap states that are expected 
to yield magnetic moments when interactions are turned 
on. Structures with global sublattice balance can still 
present midgap states and be prone to developing local 
magnetic order, at least in two situations: Distant defects 
with Ni of opposite sign and large voids with sufficiently 
long zigzag edges. 



V. DEFECTS IN SEMICONDUCTING 
RIBBONS: INTERACTION EFFECTS 

In this section we verify whether the physical picture 
anticipated from the non-interacting model remains true 
when the on-site Hubbard repulsions are included. As 
shown in previous section, unpaired spins appear in sub- 
lattice imbalanced structures. When \Nj\ > 1, the non- 
interacting model predicts that a shell of |JVj| degener- 
ate mid-gap states is half-filled. The maximization of 
the spin is expected when Coulomb repulsions are turned 
on, in the spirit of Hund's rule. At half-filling, the exact 
ground state of the Hubbard model for a bipartite lattice 



such as that of graphene satisfies Lieb's theorem^ which 
relates the sublattice imbalance and the ground state to- 
tal spin: 25 = \Nj\. For unbalanced structures this in- 
mediately confirms the Hund's rule scenario. In the case 
of balanced structures the ground state spin must be zero, 
but this could happen with local moments, as it happens 
on the edges of infinite graphene ribbons. 

The numerical calculations of this section are done 
with a unit cell of width W, length L — iV^-^a, where 
N x is the number of carbon atoms along an armchair 
chain, and with periodic boundary conditions along the 
x direction to avoid spureous zigzag edges. We consider 
unit cells as long as 25 nm and the typical number of 
atoms in a self-consistent calculation is 1000. Impor- 
tantly, our mean field results have the same relation be- 
tween the sublattice imbalance and ground state spin 
than the exact state, as predicted by Lieb theorem. 



A. Single void with (7/0 

We first revisit the single void samples. The ground 
state of structures with single atom vacancies have one 
unpaired electron within the U — model and, accord- 
ing to Lieb's theorem, spin one half in the finite U model. 
Our mean field calculation, for U — 2eV agrees with the 
Lieb theorem. There is a spin splitting of the midgap 
state As and a smaller spin splitting <5 of the conduc- 
tion and valence band states, as shown in Fig. [6l The 
spin degeneracy is thus broken, with only one of the spin 
channels of the midgap state occupied, the other being 
empty. This results in a finite magnetization density, lo- 
calized around the vacancy, as shown in Fig. [5) Although 
the magnetization resides mostly in the majority sublat- 
tice, interactions induce some reversed magnetization in 
the other sublattice. 



1. Analytical model 

We can gain some insight by doing an analytical de- 
scription of the mean field results that involves some 
approximations valid when U is much smaller than the 
band-gap of the ideal ribbon E g . In this case, we assume 
that only the midgap state is spin polarized: 



{m)o = ^l<M«)l I 



(7) 



where | 2 is the U = square modulus of the midgap 
wave function. Notice that the normalization of <p v (i) en- 
sures that the total spin of the ground state is consistent 
with Lieb's theorem, X^( m i)o = \- The corresponding 
exchange splitting is 




As = e T-eoi = Uj2\Mi)\ 2 



[mi)o 



(8) 
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FIG. 6: (Color online) (a) The spin resolved DOS for a ribbon 
with W — 7a and U = 2eV with one vacancy. Spin j (J.) is 
plotted as a positive (negative) number as a function of energy 
(we take the Fermi energy as zero), (b) Zoom of the spin-split 
midgap state, (c) Zoom of the conduction band minima. For 
clarity, we substitute the delta functions composing the DOS 
by gaussian functions with a finite broadening. 



We see that, within the simplified analytical model, the 
spin splitting of the midgap state is proportional to the 
inverse participation ratio, 77 = J^. |^>u (*) | 4 This quan- 
tity measures the degree of localization of the zero-energy 
state. An extended state in which the wave function is 
equally shared by N atoms, has 77 = -jk?- In the oppo- 
site limit where the state is localized in a single atom we 
would have 77 = 1. The inverse participation ratio shown 
in Fig. EJc) corresponds to a number of atoms in the 
range N — 5 to N = 9. As discussed above, the localiza- 
tion of the midgap states plays an important role in the 
minimal distance at which they are effectively decoupled. 

In Fig. E£a) we plot the U = gap of the ideal ribbon 
Eg and the U = 2eV spin splitting A5 of the mid-gap 
state as a function of the ribbon width W, as obtained 
from the full numerical calculation. As discussed above, 
we exclude the widths that give E g = 0. We see that 
the midgap spin splitting is a decreasing function of W. 
This is related to the fact that, in the small U limit, As is 
proportional to the inverse participation ratio ?y, which is 
also a decreasing function of the ribbon width, as shown 
in Fig. [T^c). The extension of the midgap state increases 
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FIG. 7: (Color online) (a) Non-interacting (U = 0) energy 
gap (circles) and interacting (U — 2eV) midgap spin splitting 
(squares) as a function of the ribbon width W. (b) Midgap 
spin splitting as a function of U for two ribbon widths W = 7a 
and W — 13a. (c) U = inverse participation ratio for the 
midgap state, (d) Standar deviation of the magnetization £, 
as defined in Eq. [9] as a function of U for two ribbon widths. 



as the ribbon becomes wider, resulting in a reduction of 
the midgap spin splitting. As E g tends to zero (bulk 
graphene) the midgap state becomes non-normalizable 51 
and Ag is expected to vanish (see below). 

Whereas the total magnetic moment X)i( m «) ^ s given 
by the sublattice imbalance, the degree of localization of 
the spin texture is not. In order to quantify it we define 
the standard deviation 



e = JEK) ; 



(9) 



pared to 0.57?, which is 0.016 for N y = 7 and 0.008 
for N y — 13. The non-linear terms arise from the 
interaction-driven mixing between the midgap states and 
the conduction states. This is also consistent with the 
fact that E increases as a function of U as shown in the 
lower panel of Fig. [7f d). The staggered magnetization 
is an increasing function of U. The coefficient of the 
cuadratic term decreases with the length of the sample, 
namely, with the distance between vacancies since we are 
using periodic boundary conditions. We will come back 
to this issue in the last section. 



2. Spin-charge separation for U > 

We have verified that the ground state of structures 
with single atom vacancies are locally neutral also with 
U > 0: The integrated electronic occupation in every site 
is one. Hence, a localized spin texture with total spin 1/2 
occurs in the absence of any charge localization. Our nu- 
merical results show that the addition of an extra electron 
to single vacancy structures results in a many electron 
state with total spin S = 0, local magnetization which 
is zero everywhere, and local charge accumulated in the 
same atoms and with the same distribution than the mag- 
netic texture of the charge neutral structure. These re- 
sults are shown in Fig. [8] for a ribbon with W = 7a and 
U = 2eV. Hence, it is apparent that the monoatomic va- 
cancy results in a multielectronic state with spin charge 
separation: The neutral ground state has a net electric 
charge q = 0, but a total spin S — 1/2 localized in a non- 
homogeneous spin texture in locally neutral atoms. The 
charged ground state has a net charge q = — 1, total spin 
S = 0, no local magnetic moments, and a charge texture 
localized at the same location than the spin texture of the 
neutral ground state. This phenomenon resembles that 
reported by Su-Schrieffer-Heeger™ in polyacetylene. 



In this definition E is not normalized as usual by N, the 
total number of atoms of the sample, since E character- 
izes a localized object. For sufficiently large simulation 
cells, doubling N would imply a decrease of a normal- 
ized E without changing the local properties of the local- 
ized magnetic texture. Notice that, within the analytical 
model valid for U « E g , we have E ~ h\/V- Hence, 
in the absence of staggered magnetization, both E and 7? 
would measure the localization of the magnetic moments. 
For instance, if (m.,) > at all sites and S — 1/2, the 
maximal E would be 0.5. However, the graphene lattice 
responds with a staggered magnetization to the presence 
of defects and E also measures the magnitude of that re- 
sponse. In Figs. E£b) and (d) we plot the mid-gap spin 
splitting As and E for two ribbons with W — 7a and 
W = 13a as a function of U. The midgap spin splittings 
can be fitted to A s {N y = 7){U) = 0.0172C/ + 0.0082C/ 2 
and to A s {N y = 13) (Z7) = 0.0062(7 + 0.00821676/7 2 . 
According to Eq. [5] the linear coefficients should be com- 



3. Larger voids 

We have also calculated the mean field magnetic struc- 
ture for sublattice imbalanced larger voids. In Fig. [2ft>) 
we show the magnetization profile for a triangular void 
with Nj = 2. For the chosen value of U = 2 eV the 
staggered magnetization is barely visible in this scale. In 
agreement with Lieb's theorem, it has a spin S = 1, made 
out of local moments localized, mostly, on the triangle 
boundaries. This object is the somehow complementary 
of the triangular graphene islands considered recently by 
two of us 3 . Figure [9] shows the ferrimagnetic spin tex- 
ture around a rhomboidal void with imbalance charge 
Nj = 3 — 3 = 0, i.e., composed of two triangular voids 
with Nj = ±3. Local moments with (rrij) ~ 0.05, three 
times smaller than those formed in the edges of infinite 
length zigzag ribbons, are formed on opposite corners of 
the void. We have verified that, for U = 2 eV, the small- 
est void of this shape which features local moments is the 



9 




0„ ®0„ 0O 
®0 ®Q 
0, OS ovs 

®0 ®Q 

® ®® ®o 
•> os ®® 

O □ OS 
j OO ®S 

OS OO 

®0 ®® 

W0 ®0 

© OS OS 

®0 ®0 
®0 0® 



30 



O OO OO 
OO OO 
O OO OOI 
OO OO 
O OO OO 
OO OO 
O OO OOr 
OO OO 
0„ OO OO 
OO OO 

a oa o© 

OO OO 
OO OO 

o oo o-e 



30 



FIG. 8: (Color online) Spin charge separation in single atom 
vacancy. Left column: neutral case. Right column: Charged 
case. Upper panels: Charge density qt — 1. Lower panels: 
\{nii)\/ £\ |{?7jj)|. The local charge and local spin are zero ev- 
erywhere for the neutral and charged cases respectively. The 
spin texture of the neutral case is identical to the charge tex- 
ture of the charged case. 



one of the figure. The rhomboidal void is similar to the 
hexagonal islands considered in Ref. 3 in the sense that 
both have 5 = and develop local moments if they are 
sufficiently large. 




FIG. 9: (Color online) Emerging ferrimagnetic order in a 
rhombohedral void with imbalance charge Ni = 3 — 3 = 
situated the middle of a ribbon with W = 10ai for U = 2eV. 
The largest magnetic moment per atom is {m,i} = 0.05 



B. Two vacancies with U ^ 

We now study the interaction between two magnetic 
defects with local sublattice imbalance Ni = ±1. The 
Lieb's theorem warrants that, when the sign of the 
sublattice imbalance is the same for the two defects, 
the total spin of the ground state is the sum of the 
spin of the individual defects. Hence, they are coupled 
ferromagneticall y^i ^ 7 . In contrast, if the two defects 
have opposite sublattice imbalance so that the global sub- 
lattice imbalance is zero, Lieb's theorem warrants that 
the total spin is zero. Our calculations show that this can 
happen in two different scenarios: The local magnetiza- 
tion might be zero everywhere or the two defects could 
be magnetized along opposite directions, i.e., could be 
coupled antiferromagnetically. When the defects are suf- 
ficiently far apart from each other their local electronic 
structure should be identical to that of single defects. 
Hence, the spin interaction between two magnetic de- 
fects can be either ferromagnetic or antiferromagnetic, 
as in the case of indirect exchange interactions (RKKY) 
between single site magnetic momen ta 48 ! 49 , but can also 
result in the annihilation of the local magnetic order, an 
scenario that goes beyond the RKKY picture. 

In Fig. [THlwe plot the normalized standard deviation of 
the two magnetic moments £2 for ribbons with W = 7a 
and W — 13a as a function of the defect separation. We 
normalize the computed £2 to the one corresponding to 
two independent single-defect magnetic textures, v2£i. 
Here £1 is the computed single vacancy standard devia- 
tion in the same ribbon. When the defects are sufficiently 
far away £2 must tend to \/2£i, i.e., the normalized £2 
must tend to 1. We consider the effect of the ribbon 
width W, interaction strength U, and sublattice imbal- 
ance upon the magnetic interactions between the two de- 
fects. In the case of W = 7a we show both monoatomic 
vacancies lying on the same sublattice (A+A, open cir- 
cles), whose ground state total spin is 5 = 1, and on dif- 
ferent sublattices (A+B, full circles), whose ground state 
spin is S = 0. The two curves for W = 7a are calcu- 
lated with U = 2eV. At large distances, the two defects 
become decoupled, as expected. At short distances, the 
behaviour of the magnetic texture is radically different 
for both A+A and A+B structures. In the former case £ 
is enhanced, indicating the localization of the magnetic 
texture in a smaller region. Since the total spin is 1, lo- 
cal moments survive even when the two defects are very 
close. As the separation between defects increases they 
become independent from each other and £2 = v / 2£i- 
When this happens, the energy gap between 5 = 1 and 
5 = should vanish. This is an example of rule number 
1. 

In contrast to the A+A case, the local magnetization 
of the A+B structure vanishes below a minimal distance 
D c . This is an important result. In other words, there 
is a maximal density of defects above which zero energy 
states hybridize and local magnetic moments vanish. The 
critical density depends on the energy scales of the prob- 
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lem, the single particle gap E g , controlled by the ribbon 
width, and the on-site repulsion U . For fixed U the de- 
coupling distance is definitely shorter for W = 7a than 
for W = 13a. Hence, the critical (linear) density becomes 
smaller, as the ribbon width increases. This is consistent 
with the fact that both the U = inverse participation 
ratio and the U ^ Si are decreasing functions of the 
ribbon W. The wider the ribbon, the larger the dereal- 
ization of the zero energy state. Hence, the hybridization 
between the midgap states associated with each vacancy 
survives at a larger distance for wider ribbons. Finally, 
in fig. [TO] we also show £2 for W = 13a and two values of 
U, 2 and 4 eV. The decoupling distance (critical density) 
decreases (increases) as a function of U. In other words, 
interactions drive the system magnetic, as expected. The 
A+B case is the simplest example that exemplifies rules 
number 2 and 3. 
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FIG. 10: (Color online) Normalized magnetization as a func- 
tion of the distance between vacancies for four cases. Ribbon 
with W = 7a, U = 2eV, vacancies in the same sublattice 
(open circles), same ribbon, holes in different sublattices (full 
circles), ribbont with W = 13a, holes in different sublattices, 
U — 2eV (full squares) and U = 4eV (stars) 

According to the two- vacancy calculation shown in the 
figure, the magnetic (low density, large inter-vacancy 
distance) and non-magnetic (high density, low interva- 
cancy distance) phases are separated by a crossover re- 
gion. If we take as an estimate of the critical distance 
below which local moments are quenched the distance for 
which E 2 (D c )/\/2Ei = 0.5, we find that D c = 40l for 
W = 7a and D c = 65l for W = 7a, both for U = 2eV. 
The corresponding critical linear densities n c = jj- are 
n lc (W = 7a) = 2.5 x lOW" 1 and n lc (W = 13a) = 



1.5xl0 6 cm 1 , respectively. The corresponding areal den- 
sities, n 2c = W xd are n ^c{W — 7a) = 1.4 x 10 13 cm -2 
and n2c(W — 13a) = 4.8 x 10 12 cm -2 respectively. These 
numbers should be taken as order of magnitude estimates 
of the real critical density. 

In the case of A+B pairs, the crossover from the locally 
magnetic to the non-magnetic state is similar to the one 
described in compensated graphene nanoislands^: Small 
islands are non-magnetic and larger islands have mag- 
netic edges. The critical density depends on the extension 
of the magnetization, which in turn, depends on the rib- 
bon single particle gap E g (which controls the extension 
of the U = midgap states) and on the on-site repul- 
sion U. The quenching of the local moments in the A+B 
structures is definitely related to the hybridization of the 
mid-gap states described in the non-interacting model. 
This phenomenon has an analog in zigzag ribbons. The 
midgap bands are linear combinations of top and bottom 
edge states. The hybridization is negligible in the Brilloin 
zone boundary , and is much larger in the middle. As a 
result, the exchange interaction strongly renormalize the 
zone-boundary states, opening a magnetic gap, but they 
barely change in the middle of the zone^i. 



Defective graphene ribbons as diluted magnetic 
semiconductors 



The physical picture that emerges from the previous 
discussion leads to an interesting conclusion: A semi- 
conducting graphene ribbon with a density of vacan- 
cies that induce magnetism will behave like a diluted 
(para-)magnetic semiconductor (DMS)^ 2 - provided that 
the density of defects is smaller than the critical density 
defined above (this would be an example of rule number 
3). Charged excitations will present a gap and spin ex- 
citations will not. The long ran ge f erromagnetic order 
found by Pisani et al. (see Ref. |47I ) only occurs when 
the vacancies are all in the same sublattice. It remains 
an open issue whether or not such a sublattice imbalance 
might occur in reality. Unless this can be shown, one 
should not expect long range ferromagnetic order in real 
samples. 

Interestingly, the conduction and valence bands de- 
pend on the magnetic order of the local moments, which 
might be induced by application of an external magnetic 
field. In the DMS case, the conduction and valence bands 
are exchanged coupled to the local moments, provided by 
Mn atoms. At zero field the Mn spins are randomly ori- 
ented and the average spin splitting of the bands is zero. 
Application of an external field orders the Mn spins, re- 
sulting in a finite average exchange induced spin splitting 
of the bands which is much larger than the standard Zee- 
man splitting. This is known as giant Zeeman splitting. 

The same scenario might occur in semiconducting 
graphene ribbons with magnetic vacancies. When the 
inter-defect distance is larger than the critical spacing, 
D c , the local moments are independent. Application of 
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a magnetic field aligns them and induces an exchange- 
induced splitting of the conduction and valence bands 
much larger than the intrinsic Zeeman splitting—. In the 
case of graphene ribbons with vacancies we have com- 
puted the spin splitting of the bottom of the conduction 
band 



5=Eh-E\ 



(10) 



where E a is the first level above the mid-gap state [see 
Fig. Etc)] ■ Notice that the shift of the top of the valence 
band and the bottom of the conduction band is such that 
the gap, ignoring the midgap states, is spin-independent. 
Since we consider independent defects, the calculation is 
done with a single defect per unit cell. In Fig. [Tl]we plot 
8 for a ribbon with W = 7a. In the left panel we plot 5 
as a function of the inter-defect distance considering only 
values bigger than D c for which the defects are decoupled. 
For a U = 2eV we find that 8 ranges between 5 and 15 
meV. This splitting could be obtained with an applied 
magnetic field such that g/is-B » kT, yet, i5 >> g^sB. 
As in the case of real DMS, this giant Zeeman splitting 
scales linearly with the defect density, as shown in the 
inset of the right panel of Fig. QT] Since there is a max- 
imal density above which the local moments are coupled 
and eventually they vanish, the splitting 8 can not be in- 
creased indefinitely. This phenomenon also has an analog 
in DMS: Direct antiferromagnetic coupling between Mn 
spins eventually blocks the paramagnetic coupling to the 
external field. 
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FIG. 11: (Color online) Left panel: Bottom of the conduction 
band spin splitting, 5, as a function of the defect distance and 
as a function of the linear defect density (inset) for a ribbon 
with W — 7a and U —2eV. Right panel: 8 for the same ribbon 
for a fixed defect density as a function of U. 

In the right panel of Fig. [Tl] we plot 8 as a function 
of U. For small U we find that 8 is almost linear with 
U . This can be understood in the framework of the an- 
alytical model discussed above. If we consider that, to 
lowest order, the magnetization only comes from the mid- 
gap states with wave function <j) v (i) and we compute the 



splitting of the bottom of the conduction band states to 
first order perturbation theory, we obtain: 



8 = Uj2\cj> c (i)\ 2 Y,\Mi 



(11) 



where <fr c (i) is the U = single particle state of the 
bottom of the conduction band. If we approximate 
4> c ~ l/y/~N, where N is the number of atoms in the 
unit cell, and we use the normallization condition of the 
midgap states, YJ 4 |</>„(i)| 2 = 1 then we have: 



UNy 

N 



(12) 



where N v is the number of magnetic vacancies per unit 
cell. This equation accounts also for the fact that 8 scales 
linearly with the defect density. Deviations from the lin- 
ear behaviour arise due to the magnetization that arises 
from states other than mid-gap states. 

The strong sensitivity of the conduction states of the 
defective armchair ribbon on the application of a mod- 
erate magnetic field should give rise to strong spin- 
dependent magnetotransport and magnetooptic effects, 
in analogy with DMS spintronic devices. Notice that, 
in contrast to standard Mn doped II- VI semiconductors, 
for which electrical injection of carriers results in a new 
carrier mediated couplin g 53 ! , the addition of carriers in 
this system would results in the compensation of the mid- 
gap states and the disappearance of the local moments, 
as shown in Fig. [8] 



VI. DEFECTS IN BULK GRAPHENE: 
VACANCIES 

So far we have considered the electronic structure 
of semiconducting graphene ribbons with vacancies and 
voids. As shown in Fig. [TO] the critical distance for the 
quenching of the magnetic moments increases with the 
ribbon width. An important question is whether or not 
this critical distance converges to a finite value in the 
two-dimensional limit. We have also seen in Fig. [Hd) 
how the (standard deviation of the) magnetization £ as- 
sociated with vacancies decreases as the gap of the ribbon 
decreases. In this section we address the question of what 
happens to these and other results obtained above in the 
limit of infinitely wide ribbons where the gap goes to 
zero, i.e., bulk graphene. The extrapolation to the two- 
dimensional case is not straightforward. We thus con- 
sider a new strategy. Here we consider unit cells with pe- 
riodic boundary conditions in both directions. An infinite 
graphene crystal with a unit cell formed by N y parallel 
armchair- like chains, each of them containing N x carbon 
atoms. The dimension of this unit cell is (N x ^-a, N y a) 
being a the graphene lattice parameter. We are inter- 
ested in square unit cells and therefore we consider units 
cells (N x ,N y ) of sizes (24,10), (32, 13), (40,16), (48,20), 
(60,26), and (72,31). We locate one or more vacancies 
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in the unit cell considered and, as for ribbons, we ob- 
tain the eigenvalues, eigenfunctions, and the magneti- 
zation at each place of the system, (m(i)), by solving 
self-consistently the Hamiltonian. 



A. System with vacancies of the same type 

We locate a vacancy at the center of the unit cell in 
such a way that our system describes a square lattice of 
vacancies all located on the same type of atoms, e.g., B in 
our notation. We obtain that, in agreement with Lieb's 
theorem, the ground state of the system has a magnetic 
moment 5=1/2 per vacancy. We also obtain the spin 
gap of the system As, the wavefunction <p v (i) of the first 
empty state, and the local magnetization. In the limit of 
a large unit cell As/2 and <fi v (i) should be the energy and 
the wavefunction, respectively, of an isolated vacancy. In 
two dimensions we can characterize the linear size of the 
vacancy wavefunction by its first moment 

< J R>=^|r i -r ||^(z)| 2 , (13) 

i 

where ro is the position of the vacancy and the sum is 
over all atoms of the unit cell. As discussed above, the 
bipartite character of the graphene lattice produces an 
antiferromagnetic coupling between the magnetization of 
the two sublattices of the syste m 48 i 49 . We define the 
linear size of the magnetization in each sublattice as 

Mf B) = E l^-rolK) ■ (W) 

ieA{B) 

and have opposite signs, and their sum, Md = 
+ indicates the extension of the net magnetiza- 
tion. 

In Fig. [T2] we plot the linear size of the vacancy wave- 
function and of the magnetization as function of the dis- 
tance between vacancies for different values of U. The 
first thing to note is that the size of the wavefunction 
increases linearly with the size of the unit cell, practi- 
cally independent of the value of U. This result indicates 
that the electron-electron interaction almost does not af- 
fect <p v which, when U = 0, becomes a quasilocalized 
state with weight in only one sublattice (A) and decays 
as 1/r, in agreeement with analytical results^. When 
U = only the sublattice A is magnetized and Md=M d 4 . 
However, Md is considerable smaller than (R), indicating 
that the presence of the vacancy does not just creates a 
quasilocalized state, but also modifies strongly the wave- 
functions of the states in the continuum. Notice the dif- 
ference with the ribbons where the magnetization follows 
the wavefunction for sizable confinement gaps. 

As we increase U, the magnetic texture evolves in 
such a way that their size Md, as defined in Eq. ([T3)l . 
decreases. This is accompained by an increase of the 
staggered magnetization, reflecting the antiferromagnetic 
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FIG. 12: (Color online) Linear size of the wavefunction of 
the vacancy, < R >, of the net magnetization, Md, and of 
magnetization in the sublattice A, M^, as function of the 
distance between vacancies. From top to bottom the panels 
correspond to Hubbard constant, U = 0, U = 1.5 and U — 4.5 
eV. The vacancies are located in sites type B. 



tendency of the bipartite lattice, which polarizes the sub- 
lattice B in the opposite direction that sublattice A with 
Mf > M d . This effect can be rather dramatic for mod- 
erate values of U, in Fig. [T^we show that for U = 4.5 
eV the extension of the magnetization in sublattice A is 
considerable larger than (R) and Md- Note that U = 4.5 
eV is still below the critical value V ' af ~ 5.5 eV for the 
occurrence of an antiferromagnetic instability in perfect 
grapherx o ' 55 ' 56 . Thus, a network of vacancies in the 
same sublattice would have a magnetic ground state, in 
agreement with Lieb theorem, and enhanced staggered 
magnetization, compared to perfect graphene. 

We now consider the midgap spin splitting As in two 
dimensional graphene with a finite density of vacancies in 
the same sublattice. In the previous case of semiconduct- 
ing ribbons there was a strong indication that E g > As 
for any ribbon width in the single vacancy limit. Hence, 
we might expect that As vanishes in two dimensional 
graphene. When we have a finite density of defects, As 
has also an inter-defect contribution arising from the hop- 
ping term. This mechansim is possible only if the midgap 
states have weight on the two sublattices. Midgap states 
associated with defects in the same sublattice can only be 
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FIG. 13: (Color online) (a) Midgap spin splitting as func- 
tion of the distance between vacancies and (b) weight of the 
quasilocalized wavefunction on sites B. The vacancies are lo- 
cated in the sublattice B and U = 4.5 eV. 



coupled through interaction driven sublattice mixing. As 
the gap is a product of the coupling between the magne- 
tizations induced by the Hubbard interaction, the value 
of the gap increases quadratically with U to lowest order. 
Our calculations in Fig. ITBT a) show that for U ^ 0, the 
midgap spin splitting As goes to zero as the density of 
vacancies goes to zero. In Fig. HBT b) we plot the weight 
of <fr v (i) on the sublattice where the vacancy is located. 
This quantity also tends to zero as the density of defects 
decreases. Thus, in the single impurity limit the spin gap 
goes to zero and the interacting <f> v (i) lives only on one 
sublattice. 



B. Vacancies on different sublattices: 
Compensated case 

In principle, one could expect that in real graphene 
samples the number of vacancies on sublattices A and 
B are roughly equal. In this situation Lieb's theorem 
requires that the total spin of the system should be es- 
sentially zero. From the results in the case of ribbons 
one should expect an antiferromagnetic coupling between 
the magnetic moments for small concentrations of vacan- 
cies. For large concentrations the local magnetic mo- 
ments should disappear and the sample should turn non- 
magnetic. In order to study the interaction between va- 
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FIG. 14: (Color online) Value of the critical values for the 
occurrence of a magnetic texture in a system of vacancies 
on atoms of type A and B forming two square lattice inter- 
prenetrated lattices. The inset indicates the variation of the 
standard deviation of the local magnetizations as function of 
the Hubbard constant for a system of dimension (26,60). The 
error bar is an estimation of the numerical error in the calcu- 
lations. 



cancies located in different sublattices and the local mag- 
netization in a compensated system, we locate two vacan- 
cies with total sublattice imbalance equal to zero in the 
unit cell. The A and B impurities form two interpen- 
etrated square lattices in such a way that the distance 
between impurities is maximum. We have checked that 
in agreement with Lieb's theorem the 5=0 solution is the 
ground state of the system. 

We quantify the local magnetization studying the stan- 
dard deviation of (mi), 



N 



(15) 



where the sum is over all carbon atoms and N=N X x N y 
is the number of atoms in the unit cell. In the inset of 
Fig. RBI we plot a as a function of U for a unit cell of size 
(26,60), which corresponds to a density of vacancies of 
0.5 x 10 13 cm ~ 2 . We obtain that for small values of U 
the magnetization is zero everywhere and that there is a 
critical value of the Hubbard coupling, U c , for which a 
local magnetization near the vacancies appears^. This 
critical U depends on the density of vacancies, and in Fig. 
ITU we plot U c as function of the density of vacancies. We 
obtain that U c decreases with the density of vacancies 
and from our results we conclude that the in the limit 
of zero density U c tends to zero. In that limit, each va- 
cancy hosts a spin one half texture, which is decoupled 
from the others. For high density of vacancies and mod- 
erate values of U , the kinetic energy coupling between 
the vacancies is stronger than the electronic repulsion, 
and in order to minimize the energy the system makes 
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the local magnetization zero everywhere. When the den- 
sity of impurities decreases, the kinetic energy coupling 
between vacancies decreases and eventually the impuri- 
ties become uncoupled and each of them gets a total spin 
±1/2 and the system behaves as a diluted antiferromag- 
netic system^. Unless a reason for global lattice imbal- 
ance exists (unknown to us), for realistic values of U in 
the range 1 < U < 2 eV 3 we expect that in highly irra- 
diated graphene samples the local magnetization should 
be zero across the sample. 

VII. DISCUSSION AND CONCLUSIONS 

Some of the possible limitations of our approach have 
already been mentioned in the beginning. We have left 
aside the issue of the structural stability of passivated 
vacancies and voids. Away from the edges, a single 
monoatomic vacancy might not result in a local atomic 
configuration that can be described with our model. On 
the other side, the effect of a hydrogen atom atop a car- 
bon atom on graphene is very similar to the one described 
by our models. Hence, the anomalous magnetic be- 
haviour of irradiated graphite might be related to H — C 
pairs, rather than to missing atoms. It is also important 
to signal that the mean field approximation is known to 
overestimate the appearence of magnetic order and yield 
critical values of U/t smaller than those obtained with 
methods that include quantum fluctuations^. Finally, 
we have neglected both second neighbor hopping and in- 
teratomic Coulomb repulsions. Interestingly both DFT 
and mean field Hubbard model yield a very similar de- 
scription of the magnetic behaviour of graphene islands 3 
and zigzag graphene ribbons^l. This indicates that the 
couplings neglected in the simple Hubbard model have a 
small effect on the low energy electronic structure that 
dominates the physical properties. 

We now summarize the main conclusions of this work. 
In the context of our model, the main results are: 

1. The electronic properties of the defects arising from 
the removal of atoms from graphene depend dra- 
matically on an integer number, the sublattice im- 
balance (or imbalance charge) Nj = Na — Nb, 
which counts the difference in the total number 
of atoms per sublattice removed from a perfectly 
balanced graphene lattice. Nj can take values 
0,±1,±2 etc. 

2. It can be rigorously show»2£L that the single-particle 
spectrum of a structure with sublattice imbalance 
Ni has, at least, \Nj\ midgap states per spin chan- 
nel, occupied by \Nj\ electrons in neutral graphene. 

3. Repulsive Coulomb interactions will result in a 
many-body ground state with 2S — \Nj\. This is 
an exact result in the case of the Hubbard model 3 -^. 

4. Whereas the global electronic structure of a given 
graphene system is given by the Lieb, the local 



structure is not. By assigning local sublattice im- 
balance numbers to defects, provided that they are 
sufficiently apart, a set of rules to predict basic fea- 
tures of the magnetic structure has been proposed. 

5. We find that single voids with \Nj = 1| give rise to 
states with spin-charge separation, in the sense that 
a localized magnetic texture does not entail charge 
localization. The addition of a single electron to the 
system results in a many-body state with S — and 
the disappearance of the magnetic texture, which is 
substituted by a charge texture, as seen in Fig. [5] 
In this sense, the properties of these states are very 
similar to Su-Schrieffer-Heeger midgap states^. 

6. The addition rules for two voids with a givel local 
sublattice imbalance or imbalance charge present 
similarities with those of vortices, e.g., in supercon- 
ductors. When sufficiently apart, two voids with lo- 
cal imbalance +Nj and —Nj behave like two inde- 
pendent objects with local spin 2S — \Ni\. Below 
a certain distance they annihilate each other and 
the local magnetization vanishes (Fig. ITU)) . When 
two voids with the same sign are brought together, 
they result in a region with enhanced local magne- 
tization and spin 2S — \Nj\ + \Nj\, as seen in Fig. 

121 

7. In analogy with graphene nanoislands 3 , sufficiently 
large voids with N] = can still have local mag- 
netic moments. These can interpreted as if the 
large void with Ni = was the sum of two voids 
with ±N'j. An example of this is the rhomboid of 
Fig. [51 obtained from merging two triangular voids 
with Nj — ±3 back to back. 

8. Our results show that spin interactions between two 
magnetic defects of same \Ni \ can be of three types: 
Ferromagnetic, antiferromagnetic or annihilating. 
In the first case, the ground state spin is the sum 
of the spins of the magnetic defects when infinitely 
apart, in the second case the spin is the difference 
between those two. In the third case both the to- 
tal and the local spins are zero. Antiferromagnetic 
and annihilating couplings occur in lattices without 
global sublattice imbalance, whereas ferromagnetic 
coupling requires global sublattice imbalance. 

9. Our simulations show that, in balanced defec- 
tive structures, there is a maximal density of 
monoatomic vacancies that can sustain local mo- 
ments. In the case of ribbons, this critical density 
depends on the ribbon width. The critical den- 
sity also depends on U. A phase diagram for bulk 
graphene is provided in Fig 1141 

10. Depending on the density of vacancies, distributed 
randomly in the two sublattices, we distinguish 
three phases. In the very dilute limit the sys- 
tem is paramagnetic, with some common proper- 
ties with II- VI diluted magnetic semiconductors. 
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In the opposite high-density limit the local mo- 
ments are annihilated. The intermediate phase fea- 
tures antiferromagnetically coupled local moments. 
This would be a realization of the so-called diluted 
antifcrromagnet^ 

11. We predict giant Zeeman splitting in the case of 
semiconducting ribbons in the dilute limit. Up- 
pon application of an external magnetic field such 
that gusH > ksT, the magnetic moment of all the 
defects would point parallel to applied field. This 
would result in a interaction induced splitting in the 
band states, much larger than the ordinary Zeeman 
splitting, as seen in Fig. [TT] 

12. A ferromagnetic phase is not expected for defective 
graphene, unless vacancies occur predominantly in 
one of the two sublattices. Such an unbalanced 
situation would require further explanation. 

13. In the case of zero-gap graphene, we find that 



midgap states survive, even in the interacting 
case, in the very dilute limit (Fig. [T3]). Since 
ideal graphene is a semimetal, the thermodynamics 
properties of graphene might be dominated by this 
type of defects. 

Upon completion of this manuscript, a related work by 
O. V. Yazyev (aXiv:0802.1735) has been reported. 



VIII. ACKNOWLEDGEMENTS 

We acknowledge fruitful discussions with F. Guinea, 
B. Korgel, and P. Lopez-Sancho. This work has been 
financially supported by MEC-Spain under Grants Nos. 
MAT2007-65487, MAT2006-03741, and CONSOLIDER 
CSD2007-00010, and by Generalitat Valenciana under 
Grant No. ACOMP07/054. 



1 J. M. D. Coey, Solid State Science 7, 660 (2005). 

2 H. Ohldag, T. Tyliszczak, R. Holme, D. Spemann, P. Es- 
quinazi, M. Ungureanu, and T. Butz, Phys. Rev. Lett. 98, 
187204 (2007). 

3 J. Fernandez- Rossier and J. J. Palacios Phys. Rev. Lett. 
99, 177204 (2007). 

4 R. C. Longo and L. J. Gallego, Phys. Rev. B 74, 193409 

(2006) . 

5 W. Luo , S. J. Pennycook ,S. T. Pantelides, Nano Lett. 7, 
3134 (2007). 

6 P. Crespo, R. Litran, T. C. Rojas, M. Multigner, J. M. de la 
Fuente, J. C. Sanchez-Lopez, M. A. Garcia, A. Hernando, 
S. Penades, and A. Fernandez, Phys. Rev. Lett. 93, 087204 
(2004). 

7 Y. Yamamoto, T. Miura, M. Suzuki, N. Kawamura, H. 
Miyagawa, T. Nakamura, K. Kobayashi, T. Teranishi, and 

H. Hori, Phys. Rev. Lett. 93, 116801 (2004). 

8 P. O. Lehtinen, A. S. Foster, Y. Ma, A. V. Krasheninnikov, 
and R. M. Nieminen, Phys. Rev. Lett. 93, 187202 (2004) 

9 M. A. H. Vozmediano, M. P. Lopez-Sancho, T. Stauber, F. 
Guinea, Phys. Rev. B 72, 155121 (2005). 

10 H. Kumazaki and D. S. Hirashima, J. Phys. Soc. Jpn. 76, 
064713 (2007). 

11 Oleg V. Yazyev and Lothar Helm, Phys. Rev. B 75, 125408 

(2007) . 

12 T. Chanier, I. Opahle, M. Sargolzaei, R. Hayn, and M. 
Lannoo, Phys. Rev. Lett. 100, 026405 (2008). 

13 K. S. Novoselov et al, Science 306, 666 (2004). 

14 J. Scott-Bunch et al, Nano Lett. 5, 287 (2005). 

15 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. 

I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. 
Firsov, Nature 438, 197 (2005). 

16 Y Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Nature 
438, 201 (2005). 

17 C. Berger et al, Science 312, 1191 (2006). 

18 K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. 
Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, 
and A. K. Geim, Science 315, 1379 (2007). 



19 A. Geim and K. Novoselov, Nature Materials 6, 183 (2007). 

20 Zhihong Chen, Yu-Ming Lin, Michael J. Rooks, and Phae- 
don Avouris, Physica E: Low-dimensional Systems and 
Nanostructures, 40/2, 228, (2007). 

21 M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, Phys. 
Rev. Lett. 98, 206805 (2007). 

22 B. Ozyilmaz, P. Jarillo-Herrero, D. Efetov, D. A. Abanin, 
L. S. Levitov, P. Kim, Phys. Rev. Lett. 99, 166804 (2007). 

23 L. A. Ponomarenko, F. Schedin, M. I. Katsnelson, R. Yang, 
E. H. Hill, K. S. Novoselov, A. K.Geim. larXiv:0801.0160l 

24 J. Wu, W. Pisula, and K. Mullen, Chem. Rev., 107, 718 
(2007). 

25 A. L. Vazquez de Parga et al, Phys. Rev. Lett. 100, 
056807 (2008). 

26 R. Saito, M. S. Dresselhaus, and G. Dresselhaus, Physical 
Properties of Carbon Nanotubes, (Imperial College Press, 
London, 1998). 

27 P. K. Wallace, Phys. Rev. 71, 622 (1947). 

28 M. Inui, S. A. Trugman, and E. Abrahams, Phys. Rev. B 
49, 3190 (1994). 

29 Gerardo G. Naumis, Phys. Rev. B 76, 153403 (2007). 

30 V. M . Pereira, J. M. B . Lopes dos Santos, A. H. Castro 
Neto, larXiv:0712.0806l 

H. C. Longuet-Higgins, The Journal of Chemical Physics 
18, 265 (1950). 

32 Elliott H. Lieb, Phys. Rev. Lett. 62, 1201 (1989). 

33 K. Wakabayashi, M. Sigrist, and M. Fujita, J. Phys. Soc. 
Jpn. 67, 2089 (1998). 

34 Y.-W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 
97, 216803 (2006). 

35 Young-Woo Son, M. L. Cohen, and S. G. Louie, Nature 
(London) 444, 347 (2006). 

36 L. Pisani, J. A. Chan, B. Montanari, and N. M. Harrison, 
Phys. Rev. B 75, 064418 (2007V 

37 J. Fernandez-Rossier, larXiv:0710.3484l (Phys. Rev. B in 
press). 

38 K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dressel- 
haus, Phys. Rev. B 54, 17954 (1996). 



16 



L. Brey and H. A. Fertig, Phys. Rev. B 73, 235411 (2006). 
F. Munoz-Rojas, D. Jacob, J. Fernandez-Rossier, and J. J. 
Palacios, Phys. Rev. B 74, 195417 (2006). 
S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejon, Phys. 
Rev. B 66, 035412 (2002). 

M. P. Lopez-Sancho, M. C. Munoz, and L. Chico, Phys. 
Rev. B 63, 165419 (2001). 

J. Fernandez-Rossier, J. J. Palacios, L. Brey Phys. Rev. B. 

75, 205441 (2007). 

F. Munoz-Rojas, J. Fernandez-Rossier, L. Brey, and J. J. 
Palacios, Phys. Rev. B 77, 045301 (2008). 
W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. 
Lett. 42, 1698 (1979). 

O. Hod, J. E. Peralta, and G. E. Scuseria, Phys. Rev. B 

76, 233401 (2007); O. Hod, V. Barone, and G. E. Scuseria, 
Phys. Rev. B 77, 035411 (2008). 

L. Pisani, B. Montanari, and N. M. Harrison, 
larXiv:0710.0957l 

L. Brey, H. A. Fertig, and S. Das Sarma, Phys. Rev. Lett. 



99, 116802 (2007). 

49 Saeed Saremi, Phys. Rev. B 76, 184430 (2007). 

50 W. L.Wang, S. Meng, and E. Kaxiras, Nano Lett. 8, 241 
(2008); O. V. Yazyev, W. L. Wang, S. Meng, and E. Kaxi- 
ras, Nano Lett. 8, 766 (2008). 

51 Vitor M. Pereira, F. Guinea, J. M. Lopes dos Santos, N. 
M. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 96, 
036801 (2006). 

52 J . K. Furdyna, J. Appl. Phys. 64, R29 (1988). 

53 H. Boukari, P. Kossacki, M. Bertolini, D. Ferrand, J. Cib- 
ert, S. Tatarenko, A. Wasiela, J. A. Gaj, and T. Dietl, 
Phys. Rev. Lett. 88, 207204 (2002). 

54 J. Fernandez-Rossier and L. Brey, Physical Review Letters 
93, 117201 (2004). 

55 S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992). 

56 N. Peres, M. A. N. Araujo, D. Bozi, Phys. Rev. B 70, 
195122 (2004). 



